% =================================================================
% Schniepp Lab, 2018-2021
% Deconvolute Raman Sub Peaks within reasonable parameter ranges
% =================================================================
clear;

% Load data file: uncomment this line and put in the data source file name
% load('FileName.mat');

t=1;

% Setup initial fitting parameters: 
% 1. The number of sub peaks, peak positions, and peak FWHM are based on previous publications.

% 1.1 Peak positions
peakPos=[1168.0    1227.0    1240.0    1265.0    1289.0    1312.0    1331.0    1345.6    1368.6    1398.8    1420.2    1451.0...
         1539.0    1555.5    1586.0    1609.5    1641.5    1657.5    1670.5    1685.5    1698.5    1747.0];
      
% 1.2 Peak FWHM
peakWid=[23.0000   30.0000   20.0000   25.0000   17.0000   25.0000   20.0000   15.4213   30.6685   25.0333   25.3727   25.7242...
         23.2914   15.7831   18.7831   27.7831   17.7831   17.7831   17.7831   22.7831   30.7831   25.7831];
      
% 2. The mixing percentage of Lorentizan peak component
peakL  =[0.1000    0.1000    0.1000    0.1000    0.1000    0.1000    0.1000    0.0051    0.1446    0.2394    0.3914    0.1941...
         0.0100    0.2000    0.0100    0.0100    0.0100    0.0100    0.0100    0.0100    0.0100    0.0100];

% 3.1 The magnitudes of sub peaks in XX spectrum
parXX  =[1200.0    3500.0    1000.0    2540.0    840.00    1840.0    1240.0    540.00    1700.0    2000.0    2000.0    6054.0...
         826.00    2571.0    784.00    3440.0    1840.0    3340.0    10340     3840.0    3000.0    1054.0;...
         500.00    1073.0    1000.0    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00...
         500.00    473.00    22.000    919.00    919.00    919.00    919.00    919.00    919.00    919.00];

% 3.2 The magnitudes of sub peaks in ZZ spectrum
parZZ  =[1500.0    15000     2500.0    5000.0    1640.0    2340.0    1840.0    840.00    6000.0    9500.0    2400.0    7554.0...
         1426.0    3071.0    1084.0    6340.0    2340.0    2540.0    4340.0    3840.0    3000.0    1054.0;...
         500.00    1073.0    1000.0    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00...
         500.00    473.00    22.000    919.00    919.00    919.00    919.00    919.00    919.00    919.00];
     
% 3.3 The magnitudes of sub peaks in ZX spectrum
parZX  =[500.00    9000.0    6000.0    5440.0    1640.0    3340.0    2840.0    1340.0    3400.0    2300.0    2400.0    9554.0...
         1026.0    3071.0    1484.0    5640.0    1640.0    2240.0    2340.0    2340.0    2500.0    1054.0;...
         500.00    1073.0    1000.0    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00...
         500.00    473.00    22.000    919.00    919.00    919.00    919.00    919.00    919.00    919.00];
     
% 3.4 The magnitudes of sub peaks in XZ spectrum
parXZ  =[1000.0    12000     6000.0    7000.0    2340.0    4040.0    3340.0    1840.0    4500.0    2800.0    2400.0    10554 ...
         1126.0    4071.0    2084.0    6540.0    2340.0    2740.0    3540.0    2340.0    3000.0    1054.0;...
         500.00    1073.0    1000.0    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00...
         500.00    473.00    22.000    919.00    919.00    919.00    919.00    919.00    919.00    919.00];

% 3.5 The magnitudes of sub peaks in YY spectrum
parYY  =[1200.0    3500.0    1000.0    2540.0    840.00    1840.0    1240.0    540.00    1700.0    2000.0    2000.0    6054.0...
         826.00    2571.0    784.00    3440.0    1840.0    3340.0    8340      3840.0    3000.0    1054.0;...
         500.00    1073.0    1000.0    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00    919.00...
         500.00    473.00    22.000    919.00    919.00    919.00    919.00    919.00    919.00    919.00];
      
% Combine all the initial guessing parameters for our fitting process
parguess=[peakPos; peakWid; peakL; parXX; parZZ; parZX; parXZ; parYY];
      
% set up fitting lower boundaries
lbIII=[1163    1222    1235    1260    1284    1307    1326    1340    1363    1393    1415    1446;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0;...
       0       0       0       0       0       0       0       0       0       0       0       0];
   
lbI=[1534    1550    1580    1604    1639    1652    1665    1675    1694    1742;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0;...
     0       0       0       0       0       0       0       0       0       0];

lb=[lbIII,lbI];

% set up fitting upper boundaries
ubIII=[1173    1232    1245    1270    1294    1317    1336    1350    1373    1403    1425    1456;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       1       1       1       1       1       1       1       1       1       1       1       1;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
       Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf];
   
ubI=[1544    1560    1590    1614    1649    1662    1675    1685    1704    1752;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     1       1       1       1       1       1       1       1       1       1;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf;...
     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf     Inf];

ub=[ubIII,ubI];

% Initialize the variables needed for data bookkeeping
[a,b]=size(parguess);

parP=zeros(a,b,t);
parC=zeros(a,b,t);

Position=zeros(t,b);
Width=zeros(t,b);
L=zeros(t,b);

Area=zeros((a-3)/2,b,t);

% Fitting process
for i=1:t
    fprintf('Current Position is: %d\n',i);
    parP=parguess;

    parC(:,:,i)=calSubPeaks(parP(:,:,i),x,y,lb,ub,'Cal',1500000);

    % plot the result
    Plot_MultiPeaksMultiCurves(parC(:,:,i),originals);
    set(gcf,'position',[-1903          27        1187         958]);
    
    Position(i,:)=parC(1,:,i);
    Width(i,:)=parC(2,:,i);
    L(i,:)=parC(3,:,i);
    for j=1:b
        for m=1:(a-3)/2
            r=(m-1)*2+1+3;
            Area(m,j,i)=calArea(parC(r,j,i),parC(r+1,j,i),parC(2,j,i),parC(1,j,i),parC(3,j,i),0,3000);
        end
    end
end

% Calculate basic averaged values
avgPosition=mean(Position,1);
avgWidth=mean(Width,1);
avgL=mean(L,1);

% Calculate the mean of Area over the third axis (across different pages)
avgArea=mean(Area,3);

% Calculate basic standard deviation values
stdPosition=std(Position,0,1);
stdWidth=std(Width,0,1);
stdL=std(L,0,1);
stdArea=std(Area,0,3);